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ABSTRACT 

A   computer   simulation   study   of   the   effect    of  a 
thermalized   target   lattice   on  the   ranges   of   irradiating 
ions   has   been    made.       The   calculations   are  in   good  agree- 
ment  with   experimental   results. 

The   work   is   a   continuation  in  the   development   of 
a   computer   model   formulated  at   the  USNPGS    which 
takes   into    consideration  the  displacement   of   atoms   in 
the  target   lattice  as   well  as   inelastic    energy   losses  by 
the  primary  ion.       The   simulation  was    done  for   a   Xenon 
ion  striking  the    (100)    face  of  a   tungsten  target.      This 
thesis   attempts   to   establish   the  relative  importance  of 
thermal  and  inelastic    effects   in  the  n-body   model. 
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1.  Introduction 

A  history   of   the   study   of   radiation  damage   should 
be   subdivided   in   the   experimental,    analytical  and   simu- 
lation approach   to    this   problem.      Most    of   the   earliest 
theoretical   work   was   done   in   the  analytical    manner, 
but   recently   the  high   speed   computer   has   become  a 
useful   tool,    and  the   simulation  treatment   of   radiation 
damage  has   emerged  as   an  important  adjunct   to   the 
other   theoretical  approaches. 
Theoretical 

An  excellent  and  very  readable  review  of  the 
analytical  theory  is  contained  in  Dienes  &  Vineyard 
"Radiation  Effects   in   Solids"     which   should  be   consulted 

by  all   serious   students    of   radiation  damage.      In  it,    they 

p 

note  the   work  from   the  first   efforts   of   Bohr      in   1948 

on  the  interatomic   potential,    up  to  just  about   the  time 
the  first   work   was  being   done   with  a   computer. 

Seitz   and   Koehler^   produced  a   model   which   treated 
a   locally  excited   lattice  region  as   a   temperature   spike. 
More   work  on  a  particular   range   of   temperature  spikes, 
the  highly   excited   type,    was   done  by   Brinkman^»5. 


This    "temperature   spike"    model   provides    some   of   the 
prominent   theoretical    explanations   of   radiation   effects. 

Another   of    the   several   valuable  analytic    models   of 
displacement   production   was   originally   proposed  by 
Kinchin  and   Pease0.       They  assumed   that   all   collisions 
were  binary,    the   lattice  atoms   were   initially  at   rest, 
collisions   between   pairs    of    moving  atoms    were   unlikely, 
and   the   thermal    energy   of   the   lattice  is   negligible. 

At   this   time  early   workers   were   encountering 
difficulties   in   explaining   some  of   their   results.       The 
back  diffusion  of   moving  lattice  atoms   to  the  surface 
should  be  improbable  in  view   of   the  mean  free  path  of 
the  incident   and  lattice  ions.      Silsbee'    proposed  and 
Leibfried     developed  an   energy  chain  concept   which 
allowed  direct   momentum   transmission  along  close-packed 
directions   of  the  lattice.      With   this   concept,    many  of 
the  puzzling  results   could  be  explained  immediately. 

Gibson,    Goland,    Milgram   and  Vineyard^    were  able 
to   verify  this    energy   chain  concept   with   the  use  of  the 
digital   computer.      They  used  a   copper  atom   in  a   model 
of   the   copper   lattice  and  assumed  an  n-body,    i.e.    one 


in   which   the   primary   ion  is   interacting   with   several 
target   atoms   at   one   time.      With   Born-Mayer   potential 
functions    they   were  able   to   obtain   data   that   approxi- 
mated low    energy   experimental   results   quite   well. 

Erginsoy,    Vineyard  and   Englert1^    continued  work 
with   this   n-body    model   in   the  body   centered   cubic 
lattice,    and   showed  that   the  primary   created  an  inter- 
stitial at   some   distance  from   the  point   of  impact   and 
a  vacancey  at   the  original   site.      They  also   showed   the 
threshold   energy   for   this   phenomenon   was  highly   depen- 
dent  on  the   direction  of  impact. 

The  binary   collision  model  has  been  investigated 
extensively  by   Robinson  and  Oen"»^.      They  first   assumed 
a   randomized   model   of   a  solid  which   was   inadequate  to 
explain  the  strongly  penetrating  component   observed  by 
the   CHALK  RIVER    experimental   group,    and   later  a 
lattice  structure   model  for  the  target.      As   mentioned 
before,    a   range  dependance  on  the  initial   direction   was 
observed. 

Harrison-^*  ^  and  his   students  have  used  the  n-body 


treatment   to    make  further   studies    of    the   effect   of   a 
target   lattice   structure   on   the   depth   of   penetration, 
scattering   angle  and   recoil   angles.       Gay  and   Harrison*3 
compared   the   results    of   this   n-body   treatment    to   the 
binary   collision   model.       Harrison,    Leeds   and   Gay^ 
found  that   the   total   range   computed  from   the  n-body 
model   substantially   exceeded   that   obtained  from   the 
binary   collision   model. 
Experimental 

Among   the   earlier   workers   who   attempted  heavy 
ion   range   measurements    were   Schmitt   and   Sharp  -^ 
Powers  and  Whaling       and  van  Lint,    Schmitt  and 
Suffredini1'    who    measured  the  range  of   initial  ions   with 
energies   in   the  kev  ranges. 

Recently   J.A.Davies  and   co-workers  have   obtained 
a  great   deal   of   inf ormation18'^'2^*21   with  an  ingenious 
technique.      They  bombard  the  target   with   radioactive 
ion,    then  anodize  a   surface  layer,    strip  it   chemically  and 
measure  its   radioactivity.       By   repetition   of   the  anodizing 
and  stripping,    the  fraction  of  bullets   ions   remaining  after 
a   specified   range   could  be   obtained.       This   process  how- 


ever   is    limited   to    materials    which   could  be  anodized.      In 
their   work,    they   demonstrated   that    the   nuclear   reaction 
yields    were   extremely   sensitive   to    crystal   orientation  in 
each   of   the  fee,    bec   and   diamond   type   lattice. 

Lutz   and   Sizmann^S    obtained   similar   results   in   an 
experiment   with    Kr   on   Cu.       The  ion   energy   could   be 
varied   from   10    to    15    Kev  and   ranges    were   determined 
in   the    (111),    (110)    and    (100)    directions.      Again,    the 
channeling   effect   was   reported   along   certain   preferred 
directions. 


2.         Study   Objectives 

This   report   is   a   continuation  of   the   USNPGS 

program.      In  their   latest   paper  ^    Harrison,    Leeds   and 

Gay  are   careful   to   inform    the  reader   of   two   inherent 

errors    contained   in   their   assumptions ; 

The  bullet    moves   in  a   perfect   lattice,    -which   is 
undisturbed  by  thermal   displacement   of   the  atoms. 
Preliminary   calculations   indicate  that   the   thermal 
effects   are   too   large  to   neglect,    but   their  actual 
contribution  is   still  unknown.    ••••• 

All    "inelastic"   energy  losses   have  been  neglected. 
Moving  atoms   are  known  to   lose   energy  by  a 
"friction"    method  at   low   velocities   as   well  as  by 
the  more  familiar   resonance   effects   which  appear 
at  higher   energies. 

This  thesis,    in  attempting  to   establish   the  relative 

importance   of   these  two    errors,    makes   a   comparison  of 

20 
the  results  to   those  obtained  by   Davies   et  al.         in  an 

experiment   conducted  at   Chalk  River,    Canada. 


3.         Simulation  Model 

The  basic    model   is   that   originally   developed  by   the 
USNPGS   group   over   the   past   three   years.      A  brief 
discription   of   this    model   follows   but   for    exact   details, 

the   reader   is   referred   to   reference   13,    14   and   the  paper 

23 
by   Leeds 

The   simulation   model   consists   of   a   single   primary 
and  a   bcc   target   lattice.       The  primary   can  be  fired 
into   the  target  at  any  angle  and   can  be  aimed   toward 
any   desired  point;    however,    in  view    of   the   objectives 
of   this   thesis,    all  runs   were   made   with   the  primary 
striking   the  target   normally. 

The  target   lattice  initially  is   a  39  atom   structure 
with   each  atom   subject   to   a   small   random   displacement 
which   simulates   the  thermal   displacement.      As   the  ion 
proceeds   through   the  target,    the  target   is   built   on  in 
front   of   the  primary  ion.      This   rebuilding  of   the  target 
is  accomplished  by  two    "layers"   of   atoms   at   a   time. 
At   the   same  time  the  last   two   layers   are   stripped   off 
and   discarded.      The  target  lattice   can  be  rebuilt   on  all 
six  faces. 


The   heart   of    the   simulation   program    is   the   mechanism 
by   which   the   primary   ion  interacts    with   the   remaining 
atoms.       This    model   uses   a   double  iteration  procedure   to 
determine   the   forces    on   the  primary  as   it   advances.       The 
interatomic   potential  function  and  force  function  are  of 
the   exponential   type, 

F   =    exp(A  +   Bx) 
where  A  and  B   are   empirically  determined   constants   and 
x  is   the  atomic   separation. 

The  double  iteration  procedure  is  best   described  by 

Harrison   et  al.  ^ 

The  unbalanced  force....  (on  the  primary  ion) 
is  an  average  force  calculated  by  a   double 
iteration  procedure  as   follows:    (l)assume 
an  atom   at   position  1   with  velocity  lj    (2) 
calculate  the   total  force  on  the  atom   as   a 
result   of   all  the  other  atoms   in  the  lattice 
(this   means  normally  only  about   8-10   nearest 
atoms);    (3)    call  the  calculated  force,    force 
1,    and  use  the   equation  of   motion  to    move  the 
atom   to   a   temporary  position,    position  2; 
(4)    now   repeat   the  force   calculations      for 
position  2,    call  this   force  2;    (5)    go   back  to 
position  1,    and  use  the  average  of  force   1  and 
force  2   to    move   the  atom   to   a  new   position, 
position  3»      Procedure  1  through   5    constitutes 
one   "time   step".^ 

This   double  iteration  procedure  in  an  n-body   model 


8 


is    more  accurate   than   other   models   previously   tried   on 
the   digital   computer. 

The   model   used  to   describe  the  thermal   vibration 
of   the   target   atoms    (the  primary   was   not   considered 
to   carry  any  thermal    energy  but   rather   initially  placed 
in  a   fixed  position)    consisted  of   the  application  of   a 
random   displacement   to   the  atom   from   its    "perfect 
fixed  lattice"   position.      The   maximum   amplitude   of   this 
displacement   was   consistent   with   the  assumed  tempera- 
ture,   and  a   triangular  approximation  to   the  Gaussian 

distribution   was   used  to   determine  the  probability   of 

12 
displacement     .      The  random   displacement   was   computed 

separately  for   each  axis   to   give  a   proper   spherical 

displacement. 

The  approximation  to   the  Gaussian  is 

PCf>)  --  (2irf    expC-ip2) 

-  He  -  'ffe  -£  *e?J? 

~        O  -J?  5    /pi    <     CO 

This  distribution  is   multiplied  by  an  amplitude, 
dependent   on  the  assumed  target   temperature,    which   was 


an  input   parameter   to   the   program,    and  a   factor   of 
1//3    to    make   the   displacement   spherical   since   this 
correction  is   being   made   along   each   orthagonal   axis. 
Thus 


-  C4#)('  - 


•where   D   is   the   displacement   to   be  applied   and  A(T)    is 
the  thermal  amplitude. 

The   computer   supplies  a  random   number  between 
0   and  1,    this   is    equated  to   p//6  and  the   computed 
correction  then  applied  to   the  basic   lattice  position. 

The  subprogram    -which   supplies   the  above  random 
number   is   a   standard  program   available  in  the    "library" 
of   the   computer.      It   works   with  an  original  input   14 
digit   random   odd  number,    supplied  by   the   computer 
programer  for   the  first   run,    and  thereafter  by  the 
program   itself.      To   assure  the  thermal   correction  may 
be  applied  in  both   the  plus  and   minus   direction,    the 
computer   checks   the   eleventh   digit   of   the  input   para- 
meter,   if   this   digit   is   odd,    the   correction  is   subtracted, 

10 


and  If    even,    the   thermal   correction  is   added   to   the 
basic    position. 

In  the   true   physical   crystal,    each   atom's   displace- 
ment  influences    the  adjacent   atom's   thermal   displacement. 
Thus   it   is   improbable  that   two   adjacent  atoms   would  be 
displaced  to   their   maximum   amount   in  opposite   direction. 
Therefore  as   a   displacement  is   computed  for  a   particular 
lattice  site,    only   seven  tenths   of   it  is  applied  to   the 
correction  and   three  tenths   of   the  adjacent   atom's 
displacement   is   applied  to   the  correction. 

The  amplitude   of   displacement  is   of   course  depen- 
dent  on  the  assumed  temperature   of  the  target.      This 
amplitude   was   computed  in  the  manner   of  Wert  and 
Thomson  ^,      The  following  derivation  follows  this 
reference  with   slight   changes. 

Assume  a   central  atom   is   surrounded  by   six  atoms, 
one   on  each   end   of   the  three   orthagonal  axis.      Displace 
the  atom   such  that   only  the  distance  to   the  atoms   on 
one  axis   vary    (to   the  first  approximation).      Then  the 
change  in   energy  is 

11 


where  u  =   x  -   x0   and   z   Is   the  number   of   nearest 
neighbor  atoms.       The  factor   of    two    (2)    appears   because 
the   movement   of   the  atom    causes  both   the   energy   of   the 
of   the   displaced  atom  and  of   the  neighbor  atom   to   change 
by  an   equal  amount.      The  first   term   is   the   energy   of 
bond   with  the  right   neighbor  and  the  second  with  the 
left   neighbor. 

Since  the  function   E(x)    can  be   expanded  about  its 
minimum,    the  total   energy   change  for   small   displacements 
is 

M  »X-X.    ;   «   -.(%)( ff%)^ 

An  atom   bound  to  a   definite   site  by  the  potential- 
energy  law   of  the  above  equation  is  a  simple  harmonic 
oscillator.      Hence  the  force   "f"  on  each  atom  as  a 
function  of  its   displacement  is 

Thus  for  unit   displacement,*  is   the  force   on  the 
atom. 
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The   differential    equation  describing   the  atom's 
motion  is 


Vn  J^H/jt*      ~  ~«« 


■where   m   is   the   mass   of   the  atom.      A   solution  is 
and   therefore  the  frequency  is 

can  be  estimated  by  assuming   Hookes   Law   is   applied 


to  a   unit   cube  of   side  aQ. 

Y   =   Youngs  Modulus,    F   =e(,    the  force  necessary  to 

stretch  the  cube  a  unit  distance.      Since  Young's  Modulus 

11  p 

is  about   101    newton/m    ,  oi    is   approximately  25   newton/m. 

Classically,    the   energy   of   the  atom  is 

A  £    =  3  KT 

Thus  at   room   temperature 

A  E    ^      3  K* ZfO°  I 

Since 

A  E    r    [«/z )  *  E 


«  =  jssr 
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u    - 


2  *  /,  2   *   JO'2" 


2.T 

-  // 


LA       ~  Jo 

A   slightly   more   refined   calculation  gives   u  =.05 
lattice   units. 

Similarly  at   freezing  temperature  and  melting 
temperatures,    the  amplitude  is    .025   and    .08   respectively, 
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2+.         Procedure 

This    simulation   study   was   conducted  on  the   CDC 
1604    computer,    using   both   FORTRAN   60    and   symbolic 
computer   language. 

The   computer   sets   up  a  bcc   lattice   of  atoms   of 
tungsten,    see  figure   1,    and   starts   a   Xenon  ion  into   it. 
A  brief   study   of   the   symmetry   of   the  target   lattice 
shows   that   the  indicated  impact   triangle  covers  all 
possible  points   in  the    (100)    face  of   tungsten.      Within 
this   impact   triangle,    twelve  points   were  chosen, 
figure  2,    as   being   representative   of   the   triangle.       Con- 
sideration of   the  physics   of   the  problem   as   well  as   the 
economics   of   computer   time   were  both  instrumental  in 
determining  just  how   many  points   should  be  chosen  in 
this   impact   triangle. 

The  primary  ion  was   initially  fixed  one  lattice  unit 
away  from   the  target    (perpendicular   to   the  plane  of 
the  paper  in  figure  2 )    and  given  an  energy  varying 
from   1  to   100    Kev. 

Early   runs    were   made   on  position  21,    figure  2, 
since  this   the   center   of  the   channel  between  the  two 
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nearest   atoms   on   the  face   of   the   target.       Simulated 
thermal    effects   should  be   more   meaningful   on   these 
channel  and   edge  of   channel   positions   than   on   others 
where   the  primary   strikes   a   target   atom    essentially- 
head   on  and   rebounds.      Results   showed   this   to   be   true. 

After   making   several   total   range   runs,    with   an 
initial    energy   of   1    Kev  and   with  both   a   thermalized   and 
non-thermalized   lattice,    it    was   determined   that    cutting 
off   the  primary  ion  at   50    "time   steps"    (section  3) 
would  give  a   profile  from   which  an  accurate   DE/DX 
could  be   determined.       DE/DX   is   computed  by   dividing 
the  difference   of   the  primary  ion's   initial   energy  and  its 
total   energy  after   50   time  steps  by  the  depth   of   pene- 
tration along  the  direction  of  firing. 

Runs   were   made  to   determine   DE/DX   vice  the 
total   range   since  otherwise  the  computer   time  would  be 
prohibitive,    and  a   comparison  of   DE/DX   for   thermalized 
vs   non-thermalized  lattice  is  just  as   fruitful  and  accurate 
as  a   comparison  of  total  ranges   in  assessing  the   effects 
of  thermal   energy. 

One  hundred   runs   with   initial    energies   of   1   Kev 

16 


were   made  at    each   of   the   twelve   positions   and   the   mean 
DE/DX   and   deviation  from   this    mean   were   determined. 
From   this   data,    a   sample   of   twenty   five   runs    was 
determined   to   give  a   satisfactory   mean   DE/DX.       There- 
after all   samples   were   of   twenty   five   runs;    this    was 
convenient   from    the   viewpoint   of   actual   time   required 
on   the   computer  as   well  as   sufficiently  accurate 
physically. 
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5.         Results 

Figure  3    through   12    shows    comparisons    of   the   rate 
of    energy   dissipation  at   various   initial   energies   for   each 
of   the   several   impact   points.      Position  25   and  29   were 
not   reported  as   the  primary   ion,    striking   so   close   to   a 
lattice  atom,    produced  no   statistically   consistent   data 
other   than   the   fact   it   did   rebound   back   out    every   time. 

It   is   apparent   from   these  graphs   that   those 
impact   positions   on  the   edge  of  the  open   channels   of 
the   lattice  were  the  most  affected  by  the  thermalizing 
of  the  target.      Even  these,    once  the  initial   energy   was 
above  about   10    Kev,    were  not  affected  noticeably  by 
the  thermalizing   energy. 

For  these  impact   points   of   most  interest,    data 
was  also   obtained  for  lattice  temperature  of  about 
absolute   zero  and  for  the  melting  point   temperature. 
This   data   was    only  obtained  for  initial  ion   energies   of 
1   to   10    Kev,    the  most   significant   range  for   this    effect. 

The  deviation  of   each  position's   mean   DE/DX   and 
each  position's    mean  penetration  were   recorded.      Again, 
only  if  the  point   of  impact   was   close  to  an  atom   of 
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the   target   did   this   deviation  rise  appreciably   above   0 .  5# 
of   the   principal   value. 

Finally   the   total   range   of   the   5   and  20    Kev   ions 
■was   obtained  by   a   graphical   integration.       The   reciprocal 
of   the   rate   of    energy   dissipation   was   plotted  against 
the   initial   ion   energy,    and   knowing  the   total   range  at 
the  low    end   of   this   curve,    the  total   range  at   any  other 
point   could  be  found  from   the  area   under   the   curve. 

These  ranges   for   the   several   points   of  impact 
were   plotted  on  the  impact   triangle  and   contour   lines 
for   various   penetrations   were  drawn,    see  figure  13   and 
12+.      Considering  the  total  area  of   the  triangle  as   unity, 
the  area  under  any  one  of   the  contour  lines   represented 
the  percent   of   particles  not   yet   stopped  at   this   range. 
This   is   the  required   data  for   the    "Davies   plot,"   figure 
15. 

On  this   plot,    experimental   work   taken  from   the 

20 
paper  by   Davies    et  al.         is  also   plotted.      However,    the 

curves   are  artificially    matched  at   the   one  per   cent   point 

since,    quantitatively,    the   total   range   of   the  ion  from 

any  as   yet   devised   computer   program    does   not    exactly 
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match    experimental   results.      In  addition,    errors   in   the 
graphical   integration,    which   would   tend   to   all   be   towards 
the   same   side,    would  be   corrected  by   this   procedure. 
One  final   interesting   point   was   noticed  although 
not    explained,    figure   16.      In  the   course  of   taking   data, 
runs   were   made   with   various   thermal   amplitudes   varying 
from   the   value   for   temperatures    well  below   the 
theoretical  absolute   sero   to   slightly  above  the  melting 
point   temperature.      The  points,    when  plotted  versus 
the  rate  of   energy   dissipation,    produced  a   curve  con- 
sistent  with   expectations    except   between  the  value  for 
theoretical   zero   temperature  and  actual   zero  amplitude. 
Here  there  was   a   sudden  dip  for  values   of  thermal 
amplitude  between  10        and  10~    .      The  value  then  rose 
back  up  at  a   thermal  amplitude  of   zero,    i.e.    non- 
thermalized  target   lattice. 
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6.         Conclusions   and  Acknowledgements 

Thermalizing   the   target   lattice   does   appear   to 
make  a   definite   difference   in   radiation   damage   studies 
in   the   lower    energy   range.       However   this    effect   is 
slight    compared   to   the   effect   of    "inelastic"    energy- 
loses   whence   the   ion  loses   its    energy  by    "friction"   to 
the   electron  of   the   target.      In  any  future   study   of 
radiation  damage,    the   effect   of   thermal   energy   can  be 
accounted  for   by  proper   choice   of   interatomic   potential 
function  and   constants,    which  is   still   the   central  and 
most   difficult   problem    to    solve  in   the   studies   of   this 
nature. 

A   model   such  as   the  one  used  for   this   study 
appears   to   give   correct  qualitative  results   and  is 
adequate  for   comparison  to    experimental  results. 

In  future   studies   of   the  thermal   energy   effect, 
at   least   one  lesson   can  be  learned  from   this   study. 
Since   computer   running  time  is   so    expensive  and  will 
become   more   so,    any   extra   subprograms   which  are 
unnecessary   should  be  bypassed.      This   thermal   energy 
need  not   be  applied   to   the   target  for   points   of   impact 
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near   the  atoms   of   the   target   lattice,    nor   for   targets 
where   the   initial   ion   energy   is   10    Kev   or   greater. 

I    would   like   to    extend   thanks   to    Professor    D.E. 
Harrison  Jr.    for   the  help,    guidance  and  advice   offered 
during   the  preparation   of   this   paper. 
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TARGET   LATTICE    WITH   IMPACT   TRIANGLE 
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